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Bifurcations of solitary waves propagating along the interface between two ideal 
fluids are considered. The study is based on a Hamiltonian approach. It concen- 
trates on values of the density ratio close to a critical one, where the supercritical 
bifurcation changes to the subcritical one. As the solitary wave velocity approaches 
the minimum phase velocity of linear interfacial waves (the bifurcation point), the 
solitary wave solutions transform into envelope solitons. In order to describe their 
behavior and bifurcations, a generalized nonlinear Schrodinger equation describing 
the behavior of solitons and their bifurcations is derived. In comparison with the 
classical NLS equation this equation takes into account three additional nonlinear 
terms: the so-called Lifshitz term responsible for pulse steepening, a nonlocal term 
analogous to that first found by Dysthe for gravity waves and the six-wave interac- 
tion term. We study both analytically and numerically two solitary wave families 
of this equation for values of the density ratio p that are both above and below the 
critical density ratio p cr . For p > p cr , the soliton solution can be found explicitly 
at the bifurcation point. The maximum amplitude of such a soliton is proportional 
to ^/p — p cr , and at large distances the soliton amplitude decays algebraically. A 
stability analysis shows that solitons below the critical ratio are stable in the Lya- 
punov sense in the wide range of soliton parameters. Above the critical density ratio 
solitons are shown to be unstable with respect to finite perturbations. 

PACS: 05.45.Yv; 47.55.-t; 47.90.+a 



I. INTRODUCTION 

The main goal of this paper is to study bifurcations for one-dimensional internal solitary waves 
propagating along the interface between two ideal fluids with different densities p\ and pi- The 
lighter fluid with density p2 lies above the heavier fluid with density p%: p — pij p\ < 1. These 
bifurcations occur if the solitary wave velocity V coincides with the minimum phase velocity V cr 
of linear internal waves. If the upper density p2 is small enough compared to the lower density pi, 
then a bifurcation similar to that for pure gravity-capillary waves occurs [1-4]. In this case solitary 
waves undergo a supercritical bifurcation at the critical velocity: their form approaches the form of 
the envelope solitons for the focusing one-dimensional nonlinear Schrodinger equation (ID NLSE) 
[5,6]. The soliton amplitude behaves universally near the critical velocity V = V cr : it vanishes like 
(V cr — V) 1 / 2 . The width of the solitary wave increases proportionally to (V cr — V)^ 1 / 2 . 

As the density ratio p increases, the character of the nonlinear interactions changes. The four- 
wave coupling coefficient decreases and vanishes at p = p cr = (21 — 8vo)/ll [7], that is when 
p w 0.283. Such a value may not be obtained easily in a two-fluid configuration. However, it 
may be relevant in a three-layer (or more) configuration [8], or in nonlinear optics, where similar 
singularities can occur. Above the critical ratio, solitary waves undergo a subcritical bifurcation: 
at the critical velocity, their amplitude jumps from zero for p below p cr , up to finite values when 
p is above the critical density ratio. In order to describe such type of bifurcation, it is necessary 
to keep the next order terms beyond the classical ID NLSE. When the density ratio varies in a 
neighborhood of the critical density ratio, it is possible to use, as in the derivation of the classical 
NLSE, perturbation theory assuming that the interacting wave amplitudes are small. At leading 
order, one needs to keep three kinds of terms. The first one, coming from the four-wave interaction, 
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takes into account the so-called Lifshitz invariant [9]. This term is local relative to the amplitude 
of the soliton and its first spatial derivative and, as was shown in [10], appears from the expansion 
of four-wave interaction element under the assumption about its analytical dependence. However, 
for deep-water internal waves, as we show in this paper, there exists also a nonlocal term which 
has the same order of magnitude as the first one [11]. Its structure is similar to the Dysthe term 
first found for water (gravity) waves [12] (see also [13]). In the case of water waves, this term is 
responsible for the interaction of a narrow wave packet (in /c-space) with mean flow, induced by 
the packet. The third term takes into account six-wave interactions. In order to find the six-wave 
coupling coefficient, one needs to calculate all possible rcnormalizations due to three-, four- and 
five- wave interactions and therefore this partial problem requires a lot of cumbersome calculations. 
For these calculations we use the Hamiltonian formalism (see the review [14], as well as the papers 
[15,10]), which appears to be the most adequate method for this subject. In [7], a different method 
was used: the problem was reformulated as a spatial dynamical system and only the reversibility 
was exploited. It is necessary to underline also that the use of the Hamiltonian approach to study 
solitary waves gives an appropriate framework for the temporal behavior of the dynamics of solitary 
waves, e.g. their stability. Through this formalism, it is easy to perform different kinds of averaging 
and perturbations. Second, it is crucial that by applying the Hamiltonian technique the averaging 
equations of motion retain their original Hamiltonian form. In particular, this is of great help for 
the investigation of soliton stability (see, for instance, [10]). 



II. BASIC EQUATIONS 

Consider the interface z = rj(x, i) between two ideal incompressible fluids with respective densities 
pi and p2, in the presence of gravity (with the acceleration g acting down the vertical z— axis) and 
capillarity with interfacial tension a. We shall assume that the lighter fluid with density p2 occupies 
the region oo > z > 7](x, i), and respectively the heavier fluid occupies the region — oo < z < rj(x, t). 
Flows of both fluids are considered to be potential and two-dimensional. The fluid velocities are 
given by 

Vl,2 = V(/>1,2, 

where the velocity potentials 4>i and 4>2 satisfy Laplace's equation 

A^i,2 = 0. (Ill) 
These equations are subject to the following boundary conditions. Far from the interface as z — > ±oo 

<j)i,2 -> 0. 

On the interface z = f](x, t) the kinematic conditions hold: 

^ = (-v x t) x + v z ) lt2 . (II.2) 

The dynamic condition reduces to the discontinuity of pressures across the interface due to capil- 
larity: 

n. 

Pi-P2 = -°^r 




The use of Bernoulli equations in each fluid allows to rewrite the latter equation in terms of poten- 
tials and their derivatives: 
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The equations (II.l)— (II.3) conserve the total energy: 

H = K + U, (II.4) 

where the kinetic energy is equal to 



'z>jj ^ Jz<ri " 

and the potential energy is given by the expression 

U = J {pi - P2) g ~y dx + J ai^M+l-l) dx. 

As shown first in [16] (see also [17,18,14,19]), the equations of motion (II. 2) and (II. 3) together 
with the Laplace equations (II.l) represent a Hamiltonian system. The Hamiltonian coincides with 
the energy (II. 4). The new variables — (piipi — pi'fy'i) an d the interface shape r\ are canonical 
conjugate variables: 

where ^1,2 = 4>i,2\z=ri- The given Hamiltonian form generalizes Zakharov's canonical form for free- 
surface hydrodynamics [20] . A Hamiltonian formulation of the problem of a free interface between 
two ideal fluids, under rigid lid boundary conditions for the upper fluid, was also given by Benjamin 
& Bridges [21]. Craig & Groves [22] give a similar expression, by using the Dirichlet-Neumann 
operators for both the upper and lower fluid domains (see also [23]). 

The Hamiltonian can be expanded in series with respect to powers of the canonical variables. In 
this case the steepness of the interface plays the role of a small parameter of expansion. Due to the 
conservation of the total mass for each fluid this expansion begins with the quadratic term. 

It is more convenient to work in Fourier space. Let us introduce the normal variables a k by means 
of the following formulas: 




(1 + p)uj k „ 

-2jjfcp(°* -«-*)' ( IL6 ) 



1*1 

-(a k + a*_ k ), 



2(1 + pH 

where the density pi is set equal to unity and p2 = p. In these formulas 



Uk =^-±-\g(l-p)+ak 2 ]j (II.7) 

is the dispersion relation for linear internal waves and k is the wave vector directed along the x— axis 
(ID case). 

The transformation (II. 6) diagonalizes the quadratic part of the Hamiltonian, 



«,-/ 



w k \a k \ dk. 
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As a result, the equations of motion in the new variables a k take the standard form [14]: 



da k . SH 

~~dT = (IL8) 

with the Hamiltonian 

H = Ho + Hi n t, 

where H lnii is responsible for the nonlinear interactions between waves. In the given case of internal 
waves the expansion of H lnii in the wave amplitude will contain powers starting with the cubic 
terms. 

Consider a solution of Eq. (II. 8) in the form of a solitary wave propagating with constant 
velocity V. Then the potentials and the shape -q of the interface will depend on x and t through 
the combination x — Vt. In particular, the inverse Fourier transform of a k , 



V>(m) = -^f a k (ty kx dk, 



will be a function of (x — Vt) only. 

In Fourier space such dependence implies an exponential dependence in time for the normal 
variables: 

o fc (t) = c k e~ ikvt , 

where the time-independent amplitude c k is defined from the equation 

( Uk -kV)c k = -^ = f k . (11.9) 

This equation can be casted into the following variational problem: 

6(H-VP) = 0, (11.10) 

where P = J k\c k \ 2 dk is the total momentum of the wave system. This means that a solution to 
this equation represents a stationary point of the Hamiltonian H for fixed momentum P. 

A solution to Eq. (II. 9) in the form of a solitary wave is possible if the difference uj k — kV is 
sign-definite. When the equation 

UJ k = kV (11.11) 

has real roots (say k = fc ), then, in accordance with x5(x) = 0, the solution of (II. 9) will be of the 
form: 

c fc = A^ik - fc ) + /fc M/ , f ko =0. 

U) k — KV 

Hence taking the first term as the zero approximation, after iteration the solution of this equation 
can be represented as an infinite series with respect to S(k — nfc ): 

c/t = ^2 A n S(k - nk ). 

n 

In a;— space, this solution is equivalent to a periodic solution (for details, see [10,15]). Physically, 
this criterion is very transparent. The equality (11.11) is the resonance condition for Cherenkov 
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radiation of waves by an object moving with the velocity V. Due to such radiation a solitary wave 
will lose its energy and therefore cannot be steady. 

For the internal wave dispersion (II. 7) the maximum solitary wave velocity V coincides with the 
minimum phase velocity of linear waves: 

V cr = mm — . 



It occurs when 



k = fco 

At this point the values of the linear frequency and critical velocity are 



(11.12) 



lu = Lo(k Q ) = ^2Agk and V cr = ^ = y , (11.13) 



where 

A= l - p 



1 + P 

is the Atwood number. 

As the maximum solitary wave velocity is approached, the amplitude Ck given by (II. 9) reaches 
a very sharp maximum at the point k — fco, where the straight line to — kV touches the dispersion 
curve lo = uik- 



Ck 



l -J'n 2 + k Q (V cr -V) 



-l 

fk- (H.14) 



Here k = k — fco and lo" = d 2 Lo/dk 2 is the positive-definite second derivative of cok taken at fc = fco: 
lo" = wo/(2fco) > 0. Hence one can see that as V — > V cr the width of the distribution Afc tends 
to zero, which corresponds to the peak at fc = fco becoming narrower and narrower. Due to the 
nonlincaritics of the wave system this peak generates multiple harmonics near fc = nfco with integer 
n. If the amplitude of this peak is small (if, for instance, it vanishes smoothly while approaching 
the critical velocity) then we can use the perturbation theory that consists in expanding tp through 
its harmonics: 

oo 

V>(a/)= ^n{X)e mkoX \ x' = x-Vt. (11.15) 

n=— oo 

Here the small parameter 

A = y/1 - V/V cr (11.16) 

and the "slow" coordinate X = Xx' are introduced, so that ip n {X) is the amplitude of the envelope of 
n-th harmonic. The assumption that the solitary wave amplitude vanishes continuously at V = V cr 
means that the leading term of the series in Eq. (11.15) corresponds to the first harmonic, and all 
other harmonics are small with respect to the parameter A. This is the condition under which the 
nonlinear Schrodinger equation is derived (see, for example, [15,24,25]). In this case, at leading 
order in A, wc obtain the stationary NLSE (compare with [15,10,26]) 



A^'un ^^-^.i-V, - o. (ii.1T) 
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where p is related to the matrix element T koklk2k3 of four- wave interactions (see below) as 

f i = 2nT kokokoko . (11.18) 

/From now on, we drop the subscript 1 for ip\. In complete correspondence with (11.10), the 
envelope equation (11.17) can be recasted in the following variational problem: 

5{H + lu X 2 N) = 0, (11.19) 

where 



N = J \ip\ z dx 

is the number of waves or the wave action. The (averaged) Hamiltonian H is given by the expression 

H=\Jj'\^dx + H^. 
In this approximation the leading term in the interaction Hamiltonian H lrit has the form 

H W = ?Wofco J c * kC * kiCk2Ck3 S k+ki _ k2 _ k3 dhdfodhdh = | J |^| 4 dx. (11.20) 

The tilde denotes renormalization of the vertex T due to the interaction with the zeroth and second 
harmonics, corresponding to the cubic terms in the Hamiltonian H. Thus, after averaging, the 
soliton solution is a stationary point of the (averaged) Hamiltonian for fixed N. 

/From Eq. (11.17) one can see that the localized solution is possible if the coupling coefficient p 
is negative (focusing nonlinearity) . We recall that in our case w" > 0. In this case Eq. (11.17) can 
be rewritten in dimensionless variables as follows: 

-Afy + V*x + MV = 0. (11.21) 

Its soliton solution ip s is given by 

^ = — TTTT ( IL22 ) 
cosh(Ax) 

Hence it follows that while approaching the critical velocity the soliton amplitude vanishes like 
A = (1 — VjVcr) 1 ! 2 and the soliton width grows as (1 — V '/V^) -1 / 2 The latter means that our 
approximation improves when approaching the critical velocity: the wave becomes more monochro- 
matic and nonlinearity weaker. This approximation becomes exact at the critical velocity. 

As we show in the next section such a situation occurs for all interfacial solitary waves when 
the density ratio is less than the critical value p cr = (21 — 8y/E)/ll (see for example [7]). While 
increasing p the four-wave coupling coefficient p remains negative up to the critical ratio, where it 
vanishes. For p > p cr the coefficient p becomes positive, so that the nonlinear interaction in (11.17) 
changes its character, from focusing to defocusing. In this case, in order to have solitary wave 
solutions, one needs to keep next order terms beyond the classical nonlinear Schrodinger equation 
(11.17), which should provide existence of localized solutions in the form of solitary waves. Such 
solutions were computed numerically using the full water-wave equations by Laget & Dias [27]. 
Bridges ct al. [28] computed finite-amplitude travelling waves near the transition from focusing to 
defocusing. The simplest weakly nonlinear extension retains the terms due to six-wave interactions. 
Such interactions should be of the focusing type in order to compensate for the defocusing four- 
wave interaction. From this consideration it becomes clear that the soliton amplitude undergoes 
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a jump at the point V — V cr . It is easy to estimate that such a jump will be proportional to 
y//Z. In order to obtain convergence of the Hamiltonian series expansion, the jump must be small. 
In other words, such an expansion will be valid if the deviation of p from its critical value p cr is 
small enough. The appearance of the jump at the critical velocity means that the soliton undergoes 
subcritical bifurcation. Such type of bifurcation is analogous to phase transition of first order. 
If the corresponding jump is small then we have the analogue of the first-order phase transition 
close to the second-order phase transition. For phase transitions such a situation occurs in a small 
neighborhood of the so-called tri-critical point. 

Now we will give the general structure of the Hamiltonian expansion corresponding to interfacial 
waves near the critical density ratio assuming the following two dimcnsionless parameters are small: 



A = \/l — V/V cr and 9 = 1 — p/ per- 
ks mentioned before, there are in this case three main contributions to nonlinear terms (which, 
for instance, can balance dispersion, thus providing the existence of stationary solitary waves). 
Two contributions come from the expansion of the four-wave interaction Hamiltonian. Because a 
stationary localized solution is assumed to be an envelope soliton, i.e. its spectrum remains narrow 
and concentrated near k = ko, we have to expand the four- wave matrix element Tk 1 k 2 k 3 k 4 ,, keeping 
the first-order terms that are linear in Ki = hi — k . As shown in the next section, this expansion 
contains local and nonlocal terms: 

T kl k 2 k 3 k, = 7T + + «2 + «3 + «*) (H.23) 

~ K 3| + |«2 - «3| + |«2 - K4| + |«1 - «4|)- 

The constants (3 and 7 have different parity relative to reflection fco — > — fco- The coefficient (3 
changes its sign, but the coefficient 7 retains its sign under this transform. The difference in 
parities between (3 and 7 gives different contributions to the averaged four-wave Hamiltonian: 



p\^+2ii3{r x ^-^rm 2 -im 2 m 2 (n.24) 



Here k is the positive definite integral operator 

k = -d x H, 

and H is the Hilbcrt transform: 

f(x')dx' 



Hf{ X )= l -[p.V.j_ 



The Fourier transform of the kernel of the operator k is equal to |fc|. 

The third contribution, which is local in tp, corresponds to six-wave interactions: 

H {6) =~cj |Vf das, (11.25) 

where C is the corresponding coupling coefficient. 

The solitary wave shape in this case will be defined from the solution of the following variational 
problem: 

S(H + lu q X 2 N) = 0, (11.26) 
where the (averaged) Hamiltonian is given by the expression 
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H=\j u"\4, x \ 2 dx + H (i) +H (&) . (11.27) 

The terms and are denned by Eqs. (11.24) and (11.25), respectively. The variational 

problem (11.26) can be considered as resulting from averaging the problem (11.10) over 'fast' spatial 
oscillations. 

Thus, in order to solve the variational problem (11.26), we need to know four coefficients: T (= 
Tk k k k ), P, 7 and C. One can easily see that the contributions from terms proportional to /?, 7 

— (4) 

in H and the six-wave Hamiltonian can be determined independently, which makes calculations 
more simple. 



III. HAMILTONIAN EXPANSION AND MATRIX ELEMENTS 

We begin our calculations with the four-wave matrix element T1234 in order to find T and its 
"derivatives" (3 and 7. 

The usual way to calculate matrix elements consists in expanding the Hamiltonian in series with 
respect to powers of the canonical variables ^ and r\. Then one substitutes in each term ff (™> m ) 
of the Hamiltonian the variables ^ and r\ expressed in terms of the normal amplitudes and a* k 
with the help of the formulas (II. 6), and finally one symmetrizes each term ff("' m ) against all 
and at . As a result one obtains 

/n m n-\-m 

^r;.™u„ +1 ,...^ +m n °*< n < +n <^i + - + k - k n+1 ... - t n+m ) n ^. 

i=l j=l 1=1 

The Hamiltonian for interactions H mt is represented as a sum of ff (™> m ) terms: 

H int = H(n ' m) > ( IIL1 ) 

n+m>2 

and the matrix elements T^"' m ^, | fc k are symmetric with respect to all permutations inside 
of both groups of indices i = 1, n and j = n + 1, n + m. Moreover, ff("' m ) = #*("»."). 

To find the needed matrix elements, it is convenient to represent first the Hamiltonian (II. 4) in 
the following form by integration by parts: 



H 

where 



= \ J ve + (°° - p)}^ + ecr (v°°+ - °° ) dx i ( in - 2 ) 



V = Voo.e = oi ^^F - • ( IIL3 ) 



at ' 3 9§ 



t=77 



Up to the multiplier y/l + rj% , V coincides with the normal component of the velocity (vi^ ■ n) on 
the interface z — rj(x,t). The vector n = (1 + r/ x ) _1 / 2 (— r/ x , 1) is the unit normal to the interface. 

Thus, only V needs to be expressed in terms of <3/ and r\. To find this dependence we first solve 
the Laplace equations (II. 1) for <pi and fa, 

4>i,2(x, z) = exp(±zk)A lt2 (x), (HI.4) 

where Ai, 2 {x) are functions determined from the interface boundary conditions and k is the integral 
operator defined in the previous section. The operator E(zk) = exp(zfc) is defined through the 
infinite series: 
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exp(zfc) = 1 + zk + \z 2 k 2 + ^-z 3 k 3 + ■■■ . (III.5) 

The boundary values if>i t 2 of 4>i, 2 on the interface are expressed by means of the operators E(±r]k) : 

i/)i, 2 (x) = E{±f 1 k)A lt2 = exp[±r ? (a;)fc]A li2 (a;). 

Hence by calculating the derivatives of the potentials (j>i t 2(%, z) at the interface z = r/(x,t) we have 
the following expressions for Vi t 2- 



>1,2 



dz 



Vx ^hf) = { ±E ( ± Vk)-k-V x E(±rjk)^A h2 (x). 



Z=T) 

Writing down the equality between V\ and V2 yields 

|(1 + rg) E(r)k) ■ k - r, x ^E(rjk)} A 1 (x) = {- (l + rfc) E(- V k) ■ k - Vx ^-E(- V k)^ A 2 (x), 

(III.6) 

If in addition one uses the definition of 'J, 

* = E(rik)A 1 (x) - P E(-j 1 k)A 2 (x), (III.7) 

one has two relations to determine V. 
Let us introduce the operator 

G=(l + V 2 x )E(rik)-k-E(r ] k)- 1 -r ]x -^. (III.8) 

This operator represents the Green operator for one fluid (the lower fluid) that establishes the 
relation between V\ and Vi on the surface z — i](x): 



Vi = ^{l + v^E^-k-Eirjky 1 -^-^ 



ipi(x). 



With the help of G, the kinetic energy of the lower fluid K\ is expressed as follows: 

K ^ = \ J 1P1G1P1 dx. 

This relation can be taken as the definition of the Green operator G . According to this definition 
this operator is self-adjoint as it should be. Of course, this fact can be verified also by direct 
calculations, for instance, by expanding G with respect to powers of 77. In this case one needs first 
to expand the operator E^q)^ 1 : 

E^k) -1 = 1 — i]k + rjkrjk — ^fk 2 — ^^k 3 + i^vk^k 2 + \^f'k 2 r\k — rjkrjkrjk 
-^r? 4 fc 4 + KfPrfk 2 + jrikrffc + \ } r?k\k 
— -r/ 2 k 2 r]kr]k — —rfkrf l k' x r\k — -rjkrjkr] 2 k 2 + rjkrjkrjkr]k + ■■ ■ . 

Then the Green operator G is written as a series in powers of rj: 

G = G<°) + G« + G^ + G^ + G< 4 ' + • • • , (III.9) 
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where 



G (0) = k, G (1) = -ktjk - VryV, G (2) = krjkrjk - ifc (fcr? 2 + 7y 2 fc) ft, 
G (3) = ifc (fcVfcr? + 77A-77 2 fc^) fc - krjkrjkrjk - k ( ^Vr? 3 V + ^(Ar^fc J - ^k 2 rj 3 k 2 



G (4) = krjkrjkrjkrjk + ^k 2 ij 2 kr] 2 k 2 + k ^-^?7 2 fc 2 ?7 + j^ 3 *^ krjk 

+ Ifc2 C_Ifc2^4 + ^£2^ % + \% (_l v *p + n 2 PrA k 2 
4 \ 6 / 4 \ 6 / 

— ifc {ktfkr)kr) + r\kr]ktfk^ k + kijk ^—^rj 2 k 2 T] + ^rfk 2 ^ k. 

In the case of interfacial waves the Green operator G; n , defined by the relation V = <?i n O, is 
constructed by solving the linear system (111.6,111.7) by means of the operator G (III. 9) : 

G in = (G- 1 fa) + pG^i-^y 1 . (111.10) 

Here the Green operator for the upper fluid Gi is equal to —G{—r}). The latter formula means that 

G 2 = + G« - G< 2 > + G< 3 > - G< 4 > + • • • . 

Hence one can see that the Green operator G m is self-adjoint, as it should be. The total kinetic 
energy of two fluids is defined by the operator Gj n : 



K = - I i/G iM \M»-. 

The expansion of the Green operator G; n in powers of i] is expressed through the expansion of 
Gfa)- 1 : 

g far 1 = r<°> + r« + r< 2 ) + r< 3 ) + r< 4 > + • • • , 

where 

r(°) = fc- 1 , r« = -k-'G^k- 1 , r( 2 ) = fc- 1 (G( 1 )fc- 1 GW-G( 2 ))fc- 1 

r( 3 ) = fc- 1 (-G^fc^G^fc- 1 ^ 1 ) - G< 3 > + G^k-'GW + GWk-WV) k- 1 

rw = r^-G* 4 ' + g^^g^ + G^fc- 1 ^ 1 ) 

+G< 2 )fc- 1 G( 2 ) - G^^G^fe" 1 ^ 1 ) - G^k-'G^k-'G^ 
-G^k-'G^k-'G^ + G^k-'G^k-'G^k-'G^k- 1 . 

The inverse Green operator G^ 1 is the combination 

g- 1 fa) + pG-\-r]) = (p + 1) (r(°) + Ar^ + r( 2 ) + Ar< 3 > + r< 4 ' + •••). 

Hence the Green operator G m is given through the following expansion: 
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r - 4- 4- 4- 4- j. 



where 



gL 2) = (p+iy 1 (g< 2 > + (a 2 - 1^%-^) 



(in.ii) 

(111.12) 



G| n 3) = A (p + l)- 1 + (A 2 - ^G^k-'G^k-'G^ 



(III. 13) 



G\t } = (P+ l)" 1 [G (4) - ^G^fe-^ 1 ^" 1 ^ 1 ^"^ 1 ) 

+ (A 2 - 1) (G^k^G^ + G^k-'G^ - G^k-'G^k-'G^) 
This expansion of G m allows one to write down the expansion for the Hamiltonian: 

H = H + + ff W + H ^ + H^ + ■■■, 

where 



A?7 



2(P+1) 



dx, 



A 2 krjkrjk — ifc ^fcry 2 + ?? 2 fc^ fc 



*dx 



4 



g?x. 



(111.14) 



(III. 15) 



(111.16) 



(111.17) 



2(p + l) 

+ (^ 2 - 1) (d x rjd x r]k + kr]d x r\d x - d x r)kr)d x ^J 
The Hamiltonians i?( 5 ) and _ff( 6 ) are expressed respectively through G^ and G|^ (111.13,111.14): 

iJ (5) = iy tfG^tfda;, (111.18) 



40^ + / ^|,/,-. 



(III. 19) 



Note that the Hamiltonian ff( 3 ) for the interfacial case was first calculated in the paper [18]. For 
surface waves (p = or A = 1) the expression for in the form (III. 17) was presented in the 
papers [20], [29]. 

After these calculations we can find the needed matrix elements. First, let us consider the 
Hamiltonian expansion in Fourier space (III.15-III.19). For H and H^ it gives 



Ho = \j {iT^ 1 ^ \ 2 + [( l -p)9 + °k 2 ] \ Vk | 2 |dfc, 



(111.20) 
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Hi3) =-J 2^ I y[|Ai||fe 2 |+ ^1*1*2^1+2+3^. (HI.21) 

Here the subscript 1 in *i means k\ and so on, while dk 123 , = dkidk 2 dk 3 . With this notation the 
four-wave Hamiltonian takes the form 

H {4) = J Ci2|34*l*2%^l + 2+3+4^=^ - | J k 1 k 2 k 3 k 4 7] 1 r] 2 r] 3 r] 4 5 1+2+ 3 + 4 ^J^ 2 > (HI.22) 

where 

1 f A 2 

Ci2|34 = 4 ( p+ !) N l fc 2| + h\ + \h + ki\ + \k 2 + fc 3 | + | fa + fc 4 |) 

-i(A 2 - l)/cifc 2 (|fci + fc 3 | + |fci + M + \k 2 + h\ + \k 2 + k 4 \) J . 

Hence it can be checked easily that the transition to the normal variables by means of (III. 3) 
diagonalizes the quadratic Hamiltonian: 



H = J uj k \a k \ 2 dk, 

so that the equations of motion are written in the standard form (II. 8). 

Substituting now the transformation (III. 3) into (III. 21) gives the following expression for : 

# (3) = ^ J U 123 (ala* 2 a* 3 + a 1 a 2 as)S 1+2+3 dk 123 + J V^^a^ + a^a^S^-sdk^. 

Here the matrix elements U\ 23 and V1123 are: 

= wr^ I (tr£) " 2 + NN1 (IIL23) 

= { (wSlT l ~ W + " M (nL24) 

+ fe) " 2 + **" + (l£) " 1- ww + *■*■'} - 

The Hamiltonian describing 2 — > 2 interacting waves has the form 

# (2,2) = ^ y ri2|34a^a2 a 3a4^1+2-3-4*1234. 

(2 2) 

Here the primitive (non-renormalized) 4-wave matrix element 2^2134 1S §i ven by the expression 

^\/27r^ Ti2|34 = 2(i?_!_2|34 + #34|-l-2 ~ -R-13|-24 ~ -R-14|-23 (111.25) 
— -R-24|-13 — P-23|-14 + 2Pi2|34 + 2Pi3|24 + 2P23I14), 
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where 



1 ( u\U)2k,j,k/C^ 1 



a ( k 1 k 2 k 3 k 4 \ 1/2 

32(1 + ^)^ \UJlUJ2UJ3Ld4 J 

In this section we presented only three matrix elements: U123, | 2 3 and Ti 2 |34- These matrix ele- 
ments can be used not only for one-dimensional surface waves, but also after obvious generalizations 
to the general (2D) case also. All the other matrix elements can be found by the same procedure 
as the one used to find U123, Vi\ 2 3 and T 12 | 3 4. 



IV. CALCULATIONS OF COUPLING COEFFICIENTS 

However, for solitons near the bifurcation velocity V = V cr when the density ratio is close to 
the critical one, the complete knowledge of £/i2 3 , V^| 23 and T 12 | 3 4 is sufficient. Indeed we only 
need to know the coupling coefficients /i, (3, 7 and C in (11.24) and (11.25). A straightforward and 
classical way to find them is to use the diagram technique [30] based on the renormalization of 
matrix elements. We will use this approach partially, only to calculate the constant [i. The three 
other coefficients (/?, 7 and C) can be found by using the procedure of averaging with respect to 
high frequencies. In the present case, it is the carrying frequency ujq of the main harmonic. The 
amplitudes of all the other harmonics are assumed to be small in order to apply the perturbation 
technique based on the Hamiltonian expansion. 

First we compute the coefficient ^. As is well known (see [15] for example), if three- wave in- 
teractions are not resonant, they can be excluded by canonical transformations that result in the 
renormalization of the high-order matrix elements. In particular, for 2 — > 2 interacting waves, such 
a renormalization yields 

7l2|34 = 7l2|34 (IV.l) 

— ^412. 4-2^113, 1-3 ( ; 1 ; 

\W4-2+^2-^4 Wl-3 + W3 - Ldl 

— ^214,2-4^311,3-1 f ; 1" 



"^3|2,3-2^1|4, 



1-4 



CJ 2 -4 + LO4- LU 2 W3-1 + Wl — W 3 
1 1 



LJ3-2 + UJ 2 - LO3 + LO4 - U)i 

1 1 

— ^213,2-3^411,4-1 ( ; 1 ; 

v Ul 2 -3 + ^3 — ^2 W4_i + Wi - UJ4 

— ^3+4|3,4^1+2|l,2 ( h 



0->3 + 4 — LO3 — LO4 LO1+2 — U>1 — U) 2 
— U-3-4 : 3^ 4U-I-2, 1,2 ( j j h ' 



UJ3+4 + UI3 + LU4 bJi+2 +UJ1+UJ2 

Thus, in the four-wave interaction vertex, there are three contributions: the first one comes from 
G\ n ' (<~ R), the second contribution is connected with capillarity (~ P), and the last contribution 
stems from the renormalization (IV.l) due to three-wave interactions. The latter written for k\ = 

h = h = k 4 = k (T ko 

fco|fcofeo — -^0) represents the interactions with the zeroth and second 
harmonics. However the interaction with the zeroth harmonic vanishes because the three-wave 
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matrix element Vi|23 tends to zero sufficiently rapidly if one of the wavenumbers fcj tends to zero. 
Therefore 



To - T koko \ kaka + 



9V 2 on 2 

/ " v 2ka\kgka AIJ fc fc -2fco 



2bJ ka - L0 2ka 2uJ ka + L0 2ka 

From (111.23), (HI.24) and (111.25) one obtains 

7Ttt Ak Q u> ko ( k \ 1/2 

2^W 2fe0 = 2(1 + p)1/2 ^— j • 



(IV.2) 



2ttV, 



2fc |fco ko 



AkgUJ ka / kg \ 
2(l + p)V2 ^ 2fc( J 



1/2 



(V2^) 2 T, 



k Q k \k k Q 



Substituting these expressions into (IV.2) gives 

To 



5 fcp 
16(p + l)' 



2tt(1 + p)^ ^"Ir' 

where the square of the critical Atwood number A 2 cr is equal to 5/16. This gives for the critical 
value of p 



4 + V5 

in agreement with the paper [7]. For p < p cr , the four- wave coupling coefficient is negative, and 
the corresponding nonlinearity is of the focusing type. In this case, solitary waves near the critical 
velocity V cr are described by the stationary NLSE (11.17) and undergo a supercritical bifurcation 
at V = V cr [7]. For p > p cr the coupling coefficient changes sign and the bifurcation becomes 
subcritical. To find the soliton shape in this case we need to calculate three more coefficients: (3, 7 
and C . The first two coefficients are defined from the expansion of T12134 near fc» = k$. 

To find the coefficient 7 we shall use the averaging procedure. There are two contributions to 7. 
The first one comes from the Hamiltonian (III. 17). Hence one can see that nonlocality arises from 

(4) 

two terms in H : 

Nonlocal - / Y[^TV) A2 * ^ dX + f WTT) (A2 1} { ^ ' V) % ' ^ ^ (IV - 3) 

Here the brackets (• • •) denote average with respect to short-wave oscillations so that {^f}) = 
but (d^ ■ rj) ^ 0. The latter follows after substituting 77 and \& expressed in terms of the envelope 
amplitude ip: 



V = 



k 



1/2 



* = -i 



2(l + p)w 

(i + pH i1/2 



-ioj t+ik x _|__ ^* ^iuiat— ikox 



(IV.4) 



2k 



^ e -iu t+ik x _ ^p* e iuj t-ik x^ 



It is interesting to note that these approximations of the exact formulas (III. 3) have an accuracy of 
(Afc/fco) 2 , where Ak is the width of the main peak, since at the critical wavenumbcr k = k 
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dk\k ) 

In particular, computing the averages (9rj) and (89 ■ rj) from (IV. 4) gives zero for the first mean 
and 

(d*- v ) = fcolVf- 

As a result, the nonlocal Hamiltonian (IV. 3) has the form 



H 



r(4) 



(4) 

nonlocal 



i- 2 



2(p + l) 



(A 2 



'■k\ip\ 2 dx. 



(IV.5) 



It is important that i? non i oca i is negative definite. This means that this term provides focusing. 

Another contribution to 7 comes from the interaction with the zeroth harmonic, which for gravity 
waves is responsible for the mean flow induced by the wave packet. The same statement is also 
valid for interfacial waves. 

Consider the three- wave Hamiltonian (III. 16) and average with respect to fast oscillations with 
frequency iv , under the assumption that 77 and \f contain two parts: 



V 



fco 



-1 1/2 



2(1 + p)uj 

(i + pH i1/2 



^ e -iu>at+ik x _|_ ^* ^iuiot-ikox^ 



(IV.6) 



2fc 



-iioot-\-ikox 



* iujQt—ikox 



ip*e 



) +*. 



In (IV.6), rj and 4> are low-frequency quantities responsible for the mean flow induced by the high- 
frequency wave packet, concentrated at k — fco- Substituting (IV.6) into the three-wave Hamiltonian 
(III. 16) and then averaging yields 



A (3* • rf) 
(P + l) 



Afc 
+ 



\tp\ 2 V x dx. 



(IV.7) 



This Hamiltonian describes the interaction of the wave packet with the mean flow d^/dx. To 
complete the system one needs to add the quadratic Hamiltonian, 



H 



(2) 



-u/'\ip x \ 2 dx + 



^TT) +(1 ^ ) — 



where the first part relates to the dispersion of the wave packet and the second one describes long 
linear gravity waves. In this case, the total Hamiltonian is the sum of these two terms: 

Note that -Hlf (subscript LF means low frequency) is only one part of the full Hamiltonian. It is 
used to find the coefficient 7, which can be found independently. 

The Hamiltonian equations of motion for Hlf are written similarly to (II. 5) and (II. 8): 



. <5-£/lf 1 11 1 



Ak 
(P + 1) 



SH LF _ kV Ak d\j;\ 2 

~M>~ ~ (P+T) ~~ + dx ' 



LF 



"(1 - 
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Combining the last two equations gives 



1 ~ fc* Ak dlVf 



(l-p)g (p+1) (p+1) &r ' 

The left-hand side of this equation can be neglected in comparison with the right-hand side terms. 
Therefore 

ox 

This leads to the following equation for ip : 

iilH + \j'$ xx + ^fy%| 2 -1> = 0. (IV.8) 

In the context of gravity waves the nonlinear term in this equation was first found by Dysthc 
(A = 1) [12]. It is of focusing type. This corresponds to the second nonlocal contribution to the 
Hamiltonian: 

BJL«i = -/ ^-^A^kl^dx. 
Combined with (IV. 5) it gives the final form for the nonlocal interaction Hamiltonian: 

^nonlocal = ^ |#%| 2 <fo < 0. (IV.9) 

Hence the coefficient 7 is given by 

k 2 

«. _ 



This coefficient being multiplied by — 1, it corresponds to an attraction between waves (focusing 
nonlinearity) . 

Let us now consider the coefficient /?, which is responsible for the local four-wave interactions 
(11.24). It is possible to develop the same scheme as the one we used in calculating the coefficient 
7. Another way is to calculate the first derivative of the matrix element T , 12 |34 with respect to its 
arguments at the wavenumber ki = ko (of course, excluding the nonlocal interaction which has been 
defined already). We skip all these calculations and present directly the final answer: 

Ik 2 



16(1 + pY 
so that 

The six-wave coupling coefficient C is more complicated to compute. The procedure requires 
applying the diagram technique (which, in fact, is a usual perturbation theory based on the use 
of multi-scale expansion methods) . The corresponding calculations are presented in Appendix A. 
They give 

3w ' 

where M = 289 ~ 0.685961. Thus, the six-wave interaction term also leads to a focusing 
nonlinearity. 
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V. SOLITARY WAVE SOLUTIONS 



In the previous section, we computed the coefficients that are needed to analyze solitary wave 
solutions and their bifurcations near the critical density ratio. As pointed out in Section II, soli- 
tary wave solutions represent stationary points of the Hamiltonian for fixed momentum P. When 
the critical velocity is approached, the solitons are transformed into envelope solitons. The enve- 
lope solitons, like the original ones, also represent stationary points but of the mean Hamiltonian, 
averaged with respect to fast oscillations, 



H 



-I 



dx, (V.I) 



for fixed momentum. After averaging the momentum becomes the number of waves N = J \ip\ 2 dx 
multiplied by k . The variational problem (11.26) for the envelope solitons is then written as follows: 

S(H + X 2 lu N) = 0. (V.2) 

The soliton shape in this case is governed by the corresponding Lagrange-Euler equation: 

-A 2 woV + - MlVfV* + Wl 2 ^ + 7WI 2 + 3C|V|V = 0. (V.3) 

4/C 

It is important to note that all terms in (V.I) are small compared with u>oN. However, in the 
soliton solution both dispersion and nonlincarity occur at a level comparable with X 2 uj N, where 
A 2 = (V cr — V)/V cr <C 1 • It follows, for instance, from the variational problem (V.2) or from the 
equation (V.3) for the soliton shape. 

Equation (V.3) is a pseudo-differential equation. Besides differentiation it contains the integral 
operator k that introduces some complexity in its analysis. The equation (V.3) can be simplified 
by introducing the amplitude r and the phase ip of Substituting ip = re lv into Eq. (V.3) and 
separating real and imaginary parts give the following equation for the phase: 

<p x = -^r 2 . (v.4) 

Incidentally, the phase equation for the classical NLSE (11.21) is 

(p x r 2 — constant. 

Since r — > as \x\ — > 00, the constant is equal to zero and there is no dependence on x of the 
phase as opposed to the present case. In nonlinear optics, the dependence on x of the phase ip is 
called pulse chirp. By means of the relation (V.4), the phase is excluded from the equation for the 
amplitude r: 

-AW + ^%r xx - [ir 3 + 7 rfc(r 2 ) + 3Cir 5 - 0, (V.5) 

where the six-wave coupling coefficient C is renormalized as 

4-k 2 

Substituting the rescaling 

, V3w Ci , \fj,\ - A p— / 2fc 32 

* = x , r = r<— , A=— V3Ciw , 7 = 7^==== = -7= ~ 1-60603 V.6) 

2/co \n\ V 3Cl A* V3ujoG\ V397 
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in Eq. (V.5) yields 

-A 2 r + r xx - pr 3 + r 5 + ^rk(r 2 ) = 0, (V.7) 

where all primes have been omitted. The value of p is p = — 1 for p < p cr and p = 1 for p > p cr . 

From the scaling (V.6) it is seen that the new A 2 takes small values when the old A goes to 0, and, 
respectively, large values when the old \p\ approaches . Below the critical density ratio (p = —1) 
all nonlinear terms are focusing and, thus, in this case soliton solutions exist in the whole range of 
(new) A 2 . For small A 2 , the last two terms in Eq. (V.7) can be neglected and solitons are those of 
the classical NLSE: 

r= ^ A . (V.8) 
cosh(Ax) 

Thus, at the critical velocity V cr , solitons with p < undergo supercritical bifurcation: their 

1/2 

amplitude vanishes like (1 — V/V cr ) ■ When A increases, the soliton width decreases while the 
soliton amplitude grows. This behavior is confirmed with the numerical computation of solutions 
to this equation (see below). This property can also be observed with 7 = (i.e. in the absence 
of nonlocal nonlincarity) . Then this equation admits the following analytical solution expressed in 
terms of elementary functions [7], [10]: 

r 2 = 4A ' , (V.9) 

y/l + 16A 2 /3 • cosh(2Aa;) - p 

Vi + 16A 2 /3-e 2Aa; _ 



P 2 . - 
<p= -== tan 

V^l 



4A/V3 



(V.10) 



This solution is valid for both signs of p. It is seen that this solution approaches the NLS soliton 
(V.8) when A goes to zero, if p = — 1 (focusing nonlinearity) . For p = 1 the solution (V.9) at A = 
transforms into a soliton with algebraic decay at infinity: 

'* - SUITS' < v - "> 

Such solitons have been studied in [31]. 

If 7 7^ 0, the soliton solution decays exponentially as \x\ — > oo for both signs of p (V> cx e-^l), 
except for A = 0. In this case, with p = 1 (defocusing nonlinearity), an asymptotic analysis of Eq. 
(V.7) shows that the soliton solution will decay algebraically, similarly to (V.ll): 

2-,-)"^. (V.12) 

7T J \X\ 

where iVo is the number of waves on this soliton solution. This asymptotic behavior suggests to 
seek for solutions in a form similar to (V.ll): 

r 2 = -^-j, (V.13) 

where a and A are unknown constants. Substituting this ansatz into Eq. (V.7) and assuming A = 
yields 

2 3a 2 



d 2 1 A 2 ( 1 2a 2 \, rr ( 1 



(r ) = — _ 2 , .2 " , „ . .» 2 . since H 



8x v ' a yx 2 + a 2 (x^a 2 ) 2 )' \x 2 + a 2 J ~ a{x 2 + a 2 y 

A 2 A 4 



-r 2 + r 4 = - 



+ a 2 (x 2 + a 2 Y 
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(x 2 + a 2 ) 1 and (x 2 + a 2 ) 2 , one obtains two equa- 
-3a 2 + A 4 + 2a~fA 2 = 0. 



A 2 = 7 2 + 2 - 7 ^7 2 + 3. (V.14) 

by the expression 

iV = 7T (y 7 2 + 3 - 7) ~ 2.37514. (V.15) 

This is the critical soliton solution with finite amplitude at A = 0. In other words, solitons above 
the critical density ratio (/x > 0) undergo a subcritical bifurcation. The numerical value of No is 
2.38995, which is good agreement with (V.15). 

When the density ratio is equal to the critical one (ji = 0), equation (V.5) admits another 
rescaling, different from (V.6): 



x = 2Xk x, r = r 

Under this new scaling equation (V.5) becomes 

-A 2 r + r xx + -frkr 2 + r 5 = 0, (V.16) 

where tildes have been omitted and 7 is given by the last relation in (V.6). Eq. (V.16) can be 
related to the critical stationary NLS equation. When 7 = this equation is nothing else than the 
quintic NLS equation which, as is well-known (see, for instance, [32] or the review [33]), belongs 
to the class of critical NLS equations. However, this property persists in the presence of nonlocal 
nonlinearity. It follows, for instance, if one considers the behavior of the number of waves N on the 
parameter A 2 for the soliton solution. Indeed, the change 

r = V\g(Xx), £ = Ax, (V.17) 
reduces Eq. (V.16) to a form not containing A: 

-9 + .9«+7sM5 2 )+9 5 = 0. (V.18) 
Consequently, the number of waves for soliton solutions is independent on A: 

N = N cr = J g 2 (0dC 

As A — > 00, we also come up asymptotically with the same equation (V.16). In this limit, and 
independently on \x, the main balance of the term — A 2 r in equation (V.7) comes from the last two 
last terms in its left-hand side. The latter means that we have convergence of two soliton branches 
corresponding to different signs of \i as A — > 00 . 

It is easy to show also that the Hamiltonian evaluated on the soliton solution of (V.16) is equal 
to zero. However, for solitons with \x > 0, the Hamiltonian takes positive values, while for fj, < 
(focusing case) the Hamiltonian of solitons happens to be negative. This is a simple consequence 
of the variational problem (V.2). 



Setting equal terms that are proportional to 
tions for determining a and A: 

2 - A 2 - ^— = 0, 

a 

The solution of this system with a > is 

a=l (-7 + 2vV + 3), 
For this solution the number of waves is given 
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Indeed, consider the scaling transformation that keeps the number of waves N unchanged: 

r^GD' (V ' 19) 

where r is a soliton solution of equation (V.7). Under this transform the Hamiltonian H (V.l) 
becomes a function of the scaling parameter a: 



*-*/(• 



2 7^2 



1 



dx 



1 



2 3 

Because the soliton solution is a stationary point of H (V.2), 



-r dx. 



(V.20) 



dH 

da 



= 0, 



and therefore 



Hence we have 



lr 2 kr 2 



-r" ) ,/. f 



[ir^dx. 



H 1 



fir dx. 



Thus, at = 0, H s = 0. If fi < 0, £f s < and £f s > if fj, > 0. 

Solitons of both equations (V.7) and (V.16) with different values of A were found numerically 
for both signs of \i as well as for [i = 0. In order to find them, we used the method suggested by 
V.I. Petviashvili [34]. The idea of the Petviashvili method is based on the solution of diffusive type 
equations with the introduction of a new time r. Stationary solutions of this extended equation 
coincide with the soliton solution of the original equation that we are looking for. In the present 
case we developed the following algorithm. 

The iteration scheme to find solutions connects values r^ n \x) at the nth time step with A n+1 ^ 
by the relation: 



,("+!) 



= { 



r + At ■ 



-A 2 r + r xx - [ir 3 + r 5 + 7 rfc(r 2 ) M |r (rl) | . 



(V.21) 



Here At is the time step and M {r^™^} is the functional of of the form 

1/8 

M{r} = 



1 1 - A 2 r + r xx - fir 3 + -/rk(r 2 )\\ L . 2 



\r 5 \\L 2 



where ||./||l 2 is the L 2 -norm. It is seen from this equation that its stationary solution satisfies the 
SNLS equation (V.7). The presence of the norm ||r 5 ||L 2 in the denominator of M {r} provides 
attraction to the nontrivial solution of (V.7), with r ^ 0. The efficiency of this scheme has been 
demonstrated and it provides very fast convergence. The iteration scheme (V.21) works very well 
for (i > and fi = 0. For [i < 0, however, it was convenient to use in (V.21) the factor M {r} in 
the form: 



M{r} 



X 2 r 



-7rfc(r 2 )||j 



1/4 



lir 3 + r b \\ L2 
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Figures 1 and 2 show typical soliton shapes for \x — ±1 when A ^ 0. Numerically we checked that 
solitons decay at infinity like e _A ' a: '. 



o 




X 

FIG. 1. Soliton shape r(x) with parameters fi = — 1, A = 1. 




FIG. 2. Soliton shape r(x) with parameters fj, = 1, A = 1. 



For A — and \x = 1 the iteration scheme gave the soliton dependence (Fig. 3), which is in 
agreement with the analytical result (V.13,V.14) at very high accuracy (for example, the factor M 
at the stationary solution, presented in the paper, was equal to 1 within an error less than 10 -5 ). 
The time step At was equal to 0.00018. All derivatives in the equations, as well as the action 
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of the integral operator k, were calculated by means of the standard FFT program. The initial 
conditions for the iteration procedure (V.21) were taken in the form of solitons (V.9) with 7 = 0, 
which provide exponential decay while approaching the ends of the numerical interval on x. We 
chose a symmetric numerical domain [—a, a], with a = 10 for all runs except the case A = where 
we took a = 50. The soliton shape in this particular case is presented in Fig. 3. 

0.9 | 1 1 1 1 1 1 1 1 1 1 




I 1 1 1 1 1 1 1 1 1 

—25 -20 -15 -10 -5 5 10 15 20 25 

x 

FIG. 3. Soliton shape r(x) with parameters /i = 1, A = 0. 



Fig. 4 shows the dependence of N on the soliton parameter A. It is seen that the lower soliton 
branch and the upper soliton branch converge at large A. Between the upper and the lower curves 
we have a straight line corresponding to the critical solitons with N — N cr . 



o^ 



H=-1 
H=+1 

n=o 



-A — A- 



■v- -v- 
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FIG. 4. Dependence of N on A for [i = ±1, 

The same separation takes place for soliton amplitudes (Fig. 5). At the point A = solitons 
undergo bifurcations. 




FIG. 5. Dependence of soliton amplitude on A 



For negative we have a supercritical soliton bifurcation, while a subcritical occurs when \i > 0, 
with the soliton amplitude jump given by (V.14). 



VI. STABILITY OF SOLITONS 



In this section we study the stability of the solitons that we analyzed in the previous section. In 
order to do that, we use the method based on Lyapunov's theorem. Applying this theorem to the 
present Hamiltonian system means the following: a soliton considered as a stationary point of the 
Hamiltonian H for a fixed number of waves ./V is stable if it realizes a minimum (or maximum) of 
the Hamiltonian. If a stationary point represents a saddle point then the corresponding soliton is 
expected to be unstable. The latter is only an indication that solitons are unstable. For instance, 
the classical counterexample is the Hamiltonian H = p\/2 + qf/2 — p\j e l — <z|/2 for a system with 
two degrees of freedom where the saddle point p = q = is stable. 

In fact this type of indication is already available if one looks at the scaling (V.19,V.20) for which 

»<•> -(:-»?)§/•>■ 

where r s is the solitary wave solution. Hence, it is seen that for negative \i the function H (a) is 
bounded from below and its minimum, H s < 0, corresponds to the soliton. On the contrary, for 
/i > 0, this function has a maximum equal to H s > 0, which is unbounded from below as a — > 0. 
Moreover, in the latter case, it is possible to see that the stationary point corresponding to the soliton 
solution is a saddle point. It follows immediately if, in addition to the scaling transformation, one 
considers the gauge transformation ip s — > tp s e lx under which 
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/ r 2 (X,) 2 



H^H S + I rUy x ?dx. 



For /j, < 0, however, this transformation shows that the soliton solution remains the minimum 
point. We reach the same conclusion if we apply the Vakhitov-Kolokolov criterion [35] to the 
soliton solutions. The criterion states that, if 

then solitons are unstable and they are stable in the opposite case. This criterion has a simple 
physical interpretation. The quantity —A 2 is related to the energy of the soliton as a bound state. 
If by adding one particle (i.e. increasing N) this level shifts towards the continuous spectrum, 
then obviously such bound state will be unstable. As we saw in the previous section the derivative 
dN s /dX 2 is negative for \i > and becomes positive when /i < . We would like to emphasize that 
the Vakhitov-Kolokolov criterion was derived for the classical NLS equation and, strictly speaking, 
cannot be applied to our system. Thus, we have again indications that support stability of the lower 
soliton branch and, respectively, instability for the upper branch. Note that this stability indication 
is consistent with that for the cubic NLS solitons (V.8) which, as is well-known, are stable (see, e.g. 
[33]). 

Now we show that the lower soliton branch (/i < 0) is stable: solitons from this branch indeed 
realize a minimum of the Hamiltonian H for a fixed number of waves N. The dimensionless 
Hamiltonian H written in terms of amplitude and phase reads 



H = j L 2 + ^ \r 2 kr 2 ir 6 + r 2 (<p x + fir 2 ) 2 



dx. (VI. 1) 



The last term here is positive definite and vanishes exactly on the stationary solitons. Thus, we 
now should prove that the minimum of H is reached on the soliton obtained as solution of Eq. 
(V.7). 

Consider the integral 



l r 2%r 2 + V 



dx 



and get its estimate through the integrals J r x 2 dx and N. As was shown in [10], we have the 
following estimate in the absence of nonlocal terms: 

J r e dx < (JL^j J r x 2 dx, 

where N\ — J ipo 2 dx and tpo is the soliton solution for the quintic (critical) NLS equation, 

-V>o + ifoxx + 3^o = 0. 

This solution can be easily found and gives Ni — ir/2. 

For the integral Ik = J r 2 kr 2 dx, we can write the following set of inequalities: 

J r 2 kr 2 dx < max (r 2 ) J rkr dx < J (r 2 ) x dx ij r 2 dx\ ij rk 2 rdx 



< C 2 N J r x z dx. 

This inequality can be made sharper. To do that, consider the functional 
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F{r} 



J r 2 kr 2 dx 
J r 2 dx ■ J r x 2 dx 



Its minimum value will give the best constant Ci. To find it one needs to determine a minimizer 
among all stationary points of the functional F {r} . The stationary points of F {r} are defined 
from the solutions of the equation 

2r kr 2 - r + r 0xx = 0. 

The minimizer for F {r} is given by the ground soliton solution of this equation. It is a symmetric 
function without nodes. Hence the best constant 



C2,bcst — 



1 



Thus, 



I < 



7 N N 
2 TVs + 3N 



3N 2 )I'' 



dx w 1.39035. 



r 2 dx. 



(VI.2) 



This inequality allows to find the criterion for the Hamiltonian to be bounded from below. Substi- 
tuting (VI.2) into (VI. 1) for fj, < yields 



H > 



1 



7 TV 
2 7V2 



N 2 
W 2 



2 dx 



r A dx. 



(VI.3) 



Hence it is seen that the Hamiltonian is bounded from below if 

7 N N 2 

1 1 < 1 

2 N 2 + 3N 2 ~ ' 



or 



N <N 3 = 



16' 1 4 ' 7V 2 



(VI.4) 



The final step in proving the Hamiltonian boundedness is based on the estimate for the integral 
J r 4 dx. According to [15] 



/ 



r Ux < ^=N^ 2 
~ V3 



r x dx 



1/2 



Substitution of this inequality into the estimate (VI.3) results in the desired bound for H: 



H > 



N 3 
4\/3 



1 



2 No 



N 2 
3N 2 



It turns out that the numerical value of N3 — 1.3224 is almost the same as the critical number 
N cr — 1.3225, defined by solitons with /1 = 0. Thus, solitons from the lower branch which satisfy 
the criterion (VI.4) are stable not only with respect to small perturbations but also against finite 
ones. Concerning the solitons from the upper branch, they are all unstable with respect to finite 
disturbances. 

In summary, we would like to emphasize once more the difference between ./V3 and N cr although 
the derivative dN s /d\ 2 is positive for the whole lower soliton branch and according to the Vakhitov- 
Kolokolov criterion this branch is expected to be stable. Up to now it is an open question but we 
think that it is so and N 3 and N cr should coincide. 



25 



VII. CONCLUDING REMARKS 



Thus, we have analyzed the behavior of solitons near the critical density ratio p = p cr . Above 
p cr solitons undergo a subcritical bifurcation with the amplitude jump proportional to y/p — p cr . 
Therefore our theory based on the Hamiltonian expansion works when this jump is small, i.e. in a 
small vicinity of p = p cr . In order to describe solitons a generalized NLS equation is derived based 
on the Hamiltonian average over fast oscillations with the carrying frequency corresponding to the 
critical soliton velocity V = V cr . The generalized NLS equation contains three kinds of nonlinear 
terms. The first one is nothing else than the nonlinear frequency shift taking into account six-wave 
nonlinear interactions. This is the local term. Another nonlinearity is responsible for steepening 
the envelope solitons (this is the so-called Lifshitz term [9] ) . And finally we have found the nonlocal 
contribution familiar to that calculated by Dysthe for gravity waves. This Dysthe term is focusing as 
well as the six-wave nonlinear interaction term. Within the generalized NLS equation we analyzed 
the stability of solitons. In particular, we have found a region of wave intensity, N < A3, where 
solitons are stable. Their stability is based on the boundedness of the Hamiltonian from below. 
For solitons above p cr we have shown their possible instability, at least, instability with respect to 
finite perturbations. An interesting problem is the nonlinear stage of this instability, whether it 
can lead to the collapse of solitons. The latter question is important not only from the point of 
view of interfacial waves but also because there exists some similarity between the generalized NLS 
equation derived in this paper and those describing the behavior of short optical pulses in fibers 
from the femtosecond range of pulse durations. 
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APPENDIX A: SIX- WAVE COUPLING COEFFICIENT 



As pointed out in Section III, the calculation of the nonlinear coupling coefficients /3, 7 and of 
the six-wave coefficient C can be performed independently for each coefficient. Therefore for the 
coefficient C this problem reduces to the calculation of the nonlinear frequency shift for the main 
harmonic with k = fco and ujq = u)(ko) when the first contribution to the nonlinear frequency shift 
proportional to l^l 2 is equal to zero. Second, in such a case, it is enough to consider its limit of 
a monochromatic wave, instead of a quasi-monochromatic wave. This means that one needs to 
develop the perturbation theory assuming that 

a k = V2^J2C n {t)e- inuot 5(k-nk ), (A.l) 

n 

where the leading order corresponds to the main harmonic n = 1 and amplitudes of all other 
(combined) harmonics are supposed to be small in comparison with the first harmonic amplitude. 
After this remark we introduce the small parameter e so that 

d=€il>{T), (A.2) 

where T = e A t is a slow time and i/j(T) is the envelope of the main harmonic. Hence it is easy to 
see that all needed combined harmonics C n become functions of the slow time T and take in their 
expansion powers of the parameter e: 

Co = C *=e 2 Vo + eVoi, C± 2 = eV± 2 +eV±2i, CL^efy-i, C ±3 - e 3 V±3- (A.3) 

Note that in this case nonlinear interactions with the zeroth harmonic do not give any contribution 
because the corresponding matrix elements vanish by the same law as for three-wave interaction. 

The equations of motion for amplitudes C n follow from the general equation (II. 8) for normal 
amplitude a^: 

ei ~gf- + 1 M nk o) - nu) ] C n = -i g( ^f . (A.4) 

Here H mt is mean value of the interaction Hamiltonian after substitution expressions (A.l) into 
(III.l) and average. H mt contains 19 matrix coefficients needed to determine the constant C : 

Un-2, t/l2-3, Vi|_i2, Vi I — 23 7 ^-l|l-2, V2|ll, Vjj|i2, F_2|l-3, T n \ n , T n \ 2 0, 

rp rp T (l,3) ^(1,3) T (l,3) T (0,4) ^(1,4) T (2,3) T (3,3) 

J 12|12, J l-2|l-2, J l|-in> J 1|2— 21' J 3|lll' ^-3111' J 1| — 2111 ' 1 12|111' 111|111' 

where ±1 = ±fc , ±2 = ±2fc , and so on. These 19 coefficients are calculated in accordance with 
the rules explained in the third section. 

As is seen from the equation of motion (A.4) the time derivative of C n is small compared with 
the second term i [u>(nko) — ncoo] C ni except the main harmonic where this term vanishes. In its 
turn, the equation for C\ contains terms proportional to e 3 and e 5 . The terms of the third order ( ~ 
e 3 ) are nothing more than the relation (IV. 2) serving for definition of the critical Atwood number: 
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2V 2 ] n 2 ^_ 2 

11111 + o o i = °- 

2u>\ — u>2 2uj\ + LO2 

The equation of motion for -0 appears in the fifth order of e: 

^ = 2C/ H-2^>-21 + 2U 12 -,r 2 r- 3 + 2V 1 |_ 12 0-lV'2 + 2Vl h 230-2^3 

+2F_ 1 | 1 _ 2 0_lV* 2 + 2T/ 2 |110210* + 2t/ 3 | 12 3 V'2 + 2^-2,1-30-20-3 

+4T 1 ( 1 2 |2 2 >*V 2 V'o + 4T 1 ( 2 2 | ' 1 ^0*V20 + 4T 1 (2 2 2 | ) 1 _ 2 01 2 0- 2 + 3T 1 ( | 1 i 3 1 ) 11 0_^0 

+ 6T i ( i 1 2-2iV'2^-2^ + QT^i^rr-i + 6^,^02*0* 2 + srg^rr 



where all matrix elements are defined without factors of (V^r) 2 ™ ™ as they had according to the 
definition of Section III (compare with (111.23), (111.24), 111.25)). 

Here the amplitudes ipo,ipoi,ip±2,ip±2i,'4>-i,ip±3 are found with the help of equations (A. 4). 
This is a pure algebraic procedure which we performed with the help of computer. As the result of 
substitution of ip , 4>oi,' t P±2,'4 , ±2i,' t l J -i,' t P±3 into the above equation we arrive at the equation: 

where the six-wave coupling coefficient C is positive: 
The corresponding nonlinear term is focusing. 
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